Air Shower Simulations in a Hybrid Approach using Cascade Equations 
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A new hybrid approacli to air siiower simulations is described. At higliest energies, eacli particle is 
followed individually using the traditional Monte Carlo method; this initializes a system of cascade 
equations which are applicable for energies such that the shower is one-dimensional. The cascade 
equations are solved numerically down to energies at which lateral spreading becomes significant, 
then their output serves as a source function for a 3-dimensional Monte Carlo simulation of the final 
stage of the shower. This simulation procedure reproduces the natural fluctuations in the initial 
stages of the shower, gives accurate lateral distribution functions, and provides detailed information 
about all low energy particles on an event-by-event basis. It is quite efficient in computation time. 



I. INTRODUCTION 

The field of highest energy cosmic rays is an exciting 
subject with many open questions: What is the nature of 
the primary cosmic ray? What are the highest energies? 
What are possible sources/acceleration mechanisms? Is 
there clustering of events? Is there a GZK cutoff due to 
the microwave background? Ongoing (HIRES, AG AS A) 
and future (Auger, OWL, EUSO) cosmic ray experiments 
aim to shed light on these mysteries. 

At these high energies, direct measurement of the pri- 
mary cosmic ray is impossible due to the low flux, which 
is only of order one event per square-kilometer per cen- 
tury at the highest observed energies. But cosmic rays 
initiate showers in the atmosphere, a cascade of sec- 
ondary particles from collisions with air molecules, which 
themselves collide and so on. Experiments measure these 
air showers and reconstruct from their properties infor- 
mation about the primary ray at the beginning of the 
reaction. 

Air shower models are of crucial importance for the 
reconstruction of the energy and primary type. The 
straightforward approach is to model each possible in- 
teraction of hadrons, leptons and photons with air 
molecules, and trace all secondary particles. At high 
energies this leads quickly to unpractical computation 
times, since the time grows with the number of particles 
in the shower and therefore increases rapidly with the 
primary energy. A shower of even lO^^eV has more than 
10^'' particles at its maximum and would take months 
to compute. The thinning algorithm proposed by Hillas 
[1] tries to solves this problem: below a fraction /thin of 
the primary energy only a small sample of the particles 
is actually followed in detail, attributing them a higher 
weight. This procedure introduces artificial fluctuations 
and one must compromise between these and computa- 
tion time. 

People have tried to overcome these difficulties by 
defining systems of (mostly one-dimensional) transport 
equations which describe air showers [2, 3]. The numer- 
ical solutions of these equations can then be combined 
with a Monte-Carlo in order to account for natural fiuctu- 



ations due to the first interactions and for lateral spread 
of low-energy particles [4, 5]. This is the principle of the 
hybrid method. Another realization of the hybrid ap- 
proach is to use shower libraries in which presimulated 
longitudinal profiles are combined to compute the one- 
dimensional properties of air showers [6, 7]. 

In a recent paper [8], a new approach to an old idea 
was introduced: the method of cascade equations, which 
allows one to compute longitudinal characteristics of air 
showers numerically in a very short time. 

In this paper we introduce the further development 
of this approach. Traditional Monte-Carlo methods are 
combined with cascade equations in a hybrid approach. 
This allows to construct an efficient model which ac- 
counts not only for natural fiuctuations due to the first 
interactions but also for the correct 3-dimensional spread- 
ing of low-energy secondary particles. In reasonable com- 
puting time it is possible to calculate longitudinal profiles 
and lateral distribution functions with detailed knowl- 
edge about particle momenta and arrival-times, which 
are reliable on an event-by-event basis. 



II. HYBRID APPROACH TO AIR SHOWER 
MODELING 

The solution of one-dimensional cascade equations can- 
not account for natural fiuctuations or the lateral spread 
of particles. The fiuctuations can be, as already sug- 
gested, solved by doing the first interactions up to a cer- 
tain fraction / of primary energy Eo/ma. classical Monte- 
Carlo approach, where each collision is treated individu- 
ally by the chosen hadronic model. All secondary parti- 
cles below the critical energy JEq are not followed further 
on, but taken to be initial conditions for the hadronic 
cascade equations. 

The one-dimensional cascade equations are only valid 
for large enough energies that the emission angles of sec- 
ondaries can be neglected, and the whole problem can 
be treated longitudinally. Therefore the lateral spread- 
ing of particles cannot be treated in this approach and 
we return to the Monte Carlo method for the low energy 
regime. At -Emin, the output of the cascade equations 
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FIG. 1: Schematic illustration of the hybrid approach using 
cascade equations. 



- the number of particles at certain depths and ener- 
gies - is used as a source function for the Monte Carlo 
approach, by creating single particles according to the 
source function and following them individually. This 
method is able to reproduce the lateral spread of sec- 
ondary particles even though it is neglected for collisions 

with E > -Emin- 

Fig. 1 illustrates schematically the hybrid approach. 

Throughout this paper, which concentrates on estab- 
lishing the validity of the technique, we use QGSJET 
[9] as high energy hadronic model. Low energy hadrons 
are treated by GHEISHA [10]. The electromagnetic 
part is calculated by the EGS4-code system [11]. The 
bremsstrahlung and pair-production by muons is done 
using the GEANT3.21 code [12]. At this stage we ne- 
glect the LPM effect, muon-pair production and photo- 
nuclear reactions. These effects are essentially negligible 
for hadron primaries [7, 13, 14]. The program embody- 
ing this method which was used for the calculations pre- 
sented in this paper, SENECA, is available for public use 
at http: / /cosmo.nyu.edu/ ~hjdl /SENECA/. 



III. HADRONIC CASCADE EQUATIONS 

In the domain of applicability of the cascade equations, 
the shower is one-dimensional and relativistic. Therefore 
it is completely specified by where hn{E,X)dE is 
the number of particles of a given species n with energy 
in the range [E,E + dE], at an atmospheric slant depth 



X (with X = J pg_ir(x)dx ). The reaction probability of 
a particle in the atmosphere is given by its interaction 
length and decay length, so 



dK{E,X) hn{E,X) 



dX 



Xn{E) 



hn{E,X) , (1) 



where Xn{E) is the mean free path, t„ the lifetime of 

the particle and 7 its Lorentz-factor. By writing pAir = 
Po exp(— /i/Zig) = X/ho, one can rewrite these equations 
as 



dh„{E,X) hn{E,X) B„ 



K{E,X) (2) 



dX Xn{E) EX' 

where B„ is the decay constant of hadron n defined by 

Bn = mc^ho/cTn ■ (3) 

Accounting for particles produced at higher energies 
gives rise to the following system of hadronic cascade 

equations [8] 



dK{E,X) 
dX 



= -K{E,X) 



1 ^ B„ 



iXn{E) EX 

'W^n{E',E) 



(4) 



, BrnD,„,„{E',Ey 



E'X 



Xm{E') 

dE' . 



Most important are the functions Wmn{E' , E), which are 
the energy-spectra ^ of secondary particles of type n in 
a collision of hadron m with air- molecules. Dmn{E', E) 
are the corresponding decay-functions. Equation (4) is a 
typical transport equation with a source term. The first 
term with the minus-sign accounts for particles disap- 
pearing by collisions or decays, whereas the source term 
accounts for production of secondary particles by colli- 
sions or decays of particles at higher energies. The tech- 
nique for the solution is explained in detail in reference [8] 
and we discuss in the following only two major changes. 
First, the discretized functions 



' ' rr 



,{E^) = r ^Wmn{E^E')dE' (5) 



are not calculated by fitting the functions Wmn to the 
energy-spectrum given by the hadronic Monte Carlo 
model, but by direct counting of the number of parti- 
cles falling in the energy-bin defined by the limits of the 
integral in equation 5. This gives the desired precision 
as long as the number of simulated events is high, and 
avoids introducing systematic errors due to the fitting 
procedure. The binning of the discrete energies is 



Ei 



1 GeV X 10"^, 



3 



meaning n^^'^ logarithmic bins per decade. A typical 
value is 10, but as we will see below, a higher value can 
be preferable for some applications. 

Second, the equations are modified to account for an 
arbitrary atmospheric density, since a real atmospheric 
profile is somewhat more complicated than the simple 
exponential form. Since the cascade equations are solved 
in layers Xj, Xj+i = Xi + AX ,. . . with, typically, AX = 
2.5 g/cm?, one can approximate the density in each layer 
as 



PAir = ai + biX for Xi < X < Xi+i . 



(6) 



Carlo approach can be used - e.g., the EGS4 [11] pack- 
age provides a detailed Monte Carlo model for electro- 
magnetic showers in any medium but it is very time 
consuming for higher energies. We therefore apply the 
same approach as for hadronic cascades, by defining a 
system of electromagnetic cascade equations, analogous 
to (4). Due to the fact that and photons do not decay, 
the equations simplify greatly by setting the decay con- 
stants B to zero. This basically means that the showering 
is independent of altitude if one considers path lengths in 
units of g/cw?. This fact allows a further simplification: 
First one defines energy-bins by 



The parameters Oj and bi can easily be calculated from 
any function of the density. Dropping the indices, func- 
tion (2) becomes 



dhn{E,X) hn{E,X) Bn/ho 



dX Xn{E) 
which has the solution 



E{a + hX) 



hn{E,X) (7) 



hn{E,X) =Cexp(-X/A) {a + bX)' 



(8) 



where a = at and b = hi when X is in the range Xi < 
X < Xj-i-i. Defining the corresponding cascade equations 
is then straightforward. This generalization allows one 
not only to implement different atmospheres but also to 
solve for horizontal showers due to neutrino interactions 
in the atmosphere. 

The initial condition for the cascade equation for a 
particle of type m and energy Em at depth Xm is given 
by: 



hn{E, X = Xm) = Snm5{E — E„ 



(9) 



The initial condition for the cascade equation is thus in 
general a superposition of many functions like (9), which 
accounts for the natural fluctuations. 

To recapitulate, down to a certain fraction f^^'^ = 
E^^/Eq of the primary energy, all particles are followed 
with Monte Carlo method, meaning that each collision 
is simulated explicitly by the underlying event genera- 
tor. Particles falling below i?max ^^"^ filled into the initial 
condition hn{E,X). After all particles above -E'^ax ^^''^ 
processed, one can then proceed by solving the cascade 
equations. 



IV. ELECTROMAGNETIC CASCADE 
MODELING 

The electromagnetic part of the air-shower in reference 

[8] was calculated with the analytic NKG formula. This 
is advantageous for the speed of the computation but it 
has some disadvantages: one has no detailed knowledge 
of particle spectra and it introduces inaccuracies since 
the NKG formula is only an approximation. The Monte 



Ei = IGeV X 10"d,em 



the limits of each bin being 



EilO""^ <E< £;aO"d.em , 



One defines K"" 
(l^photon, 2 



as the number of particles of type n 
electron/positron) in energy bin Ej gener- 
ated in a electromagnetic shower induced by particle m of 
energy Ei , traversing a layer of air with thickness AX of 
the order of some gjcrr? . This means 1^™" = \^™"(AX) 
is a function of layer thickness AX. In our case we choose 
AX = 2.5g/cm2. 

Let g"{X) be the number of particles of type n and 
energy Ei at a given depth X. Then, 

gnX + AX)= ^ gJ^{X)Vr-^AX). 



The function VIT" 



(AX) can be calculated in reasonable 
time by the showering model EGS4 since AX is quite 
small. Once calculated, Vj^^{AX) is stored as a table, 
and can be used to calculate efficiently any electromag- 
netic shower. 



V. LOW ENERGY SOURCE-FUNCTIONS 

We wish to follow particles down to an energy Ecuti 
below which they produce a negligible signal in the de- 
tector. Air showers have a lateral expansion of secondary 
particles, so the approach of 1-D cascade equations is 
certainly wrong for calculating particles down to lowest 
energies. Therefore, wc employ the CE only to a cer- 
tain minimum energy, E^^ and E^^ for hadronic and 
electromagnetic showers respectively. These are parame- 
ters which are determined empirically as described in the 
next section. The cases of electromagnetic and hadronic 
showers are analogous, so we describe the hadronic case 
for definiteness here. 

Particles with energies E < E^'^, produced in colli- 
sions with E > E^^ contribute to the source function 
^source j-^^ ■^yj^j,^!^ jg ^j^q numbcr of produced particles 
at depth X and energy E. It obeys the equation 
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2: An example for a source function of tt . 

urcc / 
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duced particles per logarithmic energy bin dln(i5) and depth 
dX. 



%E,X) 



dx 



Wmn{E',E) 



+ 



BmDmn{E' , E) 



E'X 



dE' . (10) 



The first term from equation (4) is missing because the 
propagation of these will be described by a Monte Carlo 
method. 

An example of a typical source function of pions for 
a vertical lO^^eV proton induced shower is given in 
Fig. 2, choosing E^^ = 10^ GeV. The source func- 
tion is used to generate particles, which arc then traced 
in the Monte-Carlo part of the air shower simulation 
code. With unlimited computational speed, the num- 
ber Nn = /fSd" / — — w^—^dEdX of particles would 

cut 

be produced for each species n, however at high energies 
this is time consuming. Instead, only a certain fraction 
of the total number of particles is sampled, attributing 
to each particle a suitable weight larger than 1. A prac- 
tical way to define the sampling procedure is to spec- 
ify the total amount of hadronic (or cm) energy which 
is carried by the particles followed in the low energy 
MC. Because computation time is roughly proportional 
to energy, this procedure allows one to achieve equally 
good statistics independent of the energy Eq of the pri- 
mary cosmic ray. To be precise, the procedure is the 
following. The total energy of hadrons in the low en- 
ergy regime produced by reactions with energy greater 

than Eif^ is Ell^^,^, = Y.J^ ! E^IZiE^SldEdX, 

cut 

with the index n summing over the particle types. If 
low energy particles distributed according to the source 
function would be generated until their energy totalled 
-^lovftot' weight would be 1. Instead we produce par- 



method 


thinning level /th 


CPU-time[min] 


MC 




71 




383 


10-^ 


2148 


CE 




19 



TABLE I: CPU time comparison for the showers shown in 
Fig. 3. The showers have been calculated on a 1.266Ghz 
processor. 



tides until their total energy is E^^ < E^^^.^^, so the 
weight attributed to each particle is w = ^h^j"* > 1. 

With this method, a simple adjustment of E^^ controls 
the final weight of all particles, since the shower of each 
generated particle is followed in full detail with no further 
thinning. Thus this method also overcomes a weakness 
of the normal thinning method, where the weight itself 
can fluctuate a lot. 

It can happen that particles with E < E^^ are pro- 
duced in the high energy MC stage, in reactions with 
particles of an energy E > Ejl^^^. Following all of these 
with weight 1 down to E^^^ is time consuming and un- 
necessary for particles with angle less than about « 5° 
with respect to the shower axis. These arc stored in the 
low energy source function and re-appear in the compu- 
tation at the stage when low energy particles are created 
from the source function as discussed above. Low en- 
ergy particles with larger angles arc treated directly in 
the Monte Carlo part. 

Neutral pions have a very short decay length and are 
therefore treated separately. In the system of cascade 
equations (4), they appear only as secondary particles 
and formula (10) is evaluated for E < E^^. The result- 
ing 7r"s can then be fed either into the electromagnetic 
cascade equations for E > -E^Jfu or into the Monte Carlo 
part of the code for E < E^^. This approach is valid for 
E^f^ <B^o=3x lO^^eV. 

The generalization of the source function method to 
the electromagnetic case is straightforward. 



VI. TESTS AND APPLICATIONS 

In the ideal case all the methods described here arc just 
of technical nature, which means that the final result of 
physical observables of air showers should not depend on 
whether they are calculated with the traditional Monte 
Carlo (MC) method or with the hybrid method proposed 
here, using cascade equations (CE). Therefore a first step 
is to check the new technique by comparing the results of 
the two approaches and the influence of the parameters 
E^'^ and E^^ on longitudinal and lateral profiles. In a 
second step we show some comparisons to the CORSIKA 
model, which can be configured to use the same external 
models - QGSJET, GHEISHA and EGS4. 

Before doing so we show a comparison of the com- 
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FIG. 3: The lateral distribution function of showers calculated 
with cascade equations and Monte Carlo using thinning. The 
bottom figure shows the relative fluctuations a /p. 



putation time necessary to simulate W^^eV proton in- 
duced vertical showers in table I. The results of CE 



using E^^ 



= lO^GeV and E\ 



low 



lO^GeV are com- 



pared with a pure traditional Monte-Carlo method using 
various thinning levels /thin (/thinEo is the energy be- 
low which only one secondary i in each reaction is fol- 
lowed with probability Pi = Ei/J2i^i having a weight 

Wi = l/Pi). 

The corresponding lateral distribution functions and 

the relative fluctuations can be seen in Fig. 3. The rela- 
tive fluctuations in each lateral bin are defined by 



^/E 



(11) 



where Wi is the weight of particle i. One sees in Fig. 
3 that the quality of the LDF computed with CE is 
somewhere between the thinning levels 10~^ and 10~^ 
(approaching the latter for large distances), whereas the 
computation time is at least 20 times lower as seen in 
table I. As the energy of the primary cosmic ray is in- 
creased, the CPU time of the CE stays approximately 
constant if one keeps the same values for and Ef^. 
This is because most of the time is used for the low energy 
Monte Carlo part. In the pure MC method, higher values 
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FIG. 4: Comparison of longitudinal and lateral profiles us- 
ing the MC and CE approaches (upper and lower figures, re- 
spectively). CEl denotes cascade equations with 10 bins per 
decade, CE3 with 30. There is no noticeable difference in the 
LDF for different binnings; these are therefore not shown. 



of the thinning parameter /thin can often be used while 
maintaining the statistical quality, so the improvement 
factor using the CE grows only slowly at higher energies. 
For instance, at 3 x lO'^'^eV, /thin = 10~^ gives results 
comparable to the CE approach using the same values 
for and Ef^ with a time difference of a factor of 
about 40. 



A. Checks on an average shower basis 

If one applies the cascade equations starting from the 
primary energy, an average shower is calculated. In this 
case the initial condition consists just of the primary cos- 
mic ray. Wc can compare such a result to the average 
of many showers computed by the MC-method. Fig. 4 
shows such a comparison for lO^^eV proton induced ver- 
tical showers. The lateral and longitudinal profiles agree 
nicely within a small error. The shower maxima X^ax 
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agree within less than 1%. As for the shower size N„iax 
(number of particles at shower maximum) , we achieve 3% 
accuracy if we use 10 bins per decade in the numerical 
solution, but this can be improved to 1% by using 30 bins 
instead. 

The other relevant parameters of the CE are E"^^^ = 
lO^Gey, = lOGeV. The performance of the CE 

depends on these parameters as well as on the binning 
chosen for the numerical solution. A fine binning takes 
a long time to compute, whereas a more rough binning 
might introduce a significant error. Similarly, minimizing 
computation time argues for a low energy threshold i^min 
for both cascades, but not too low in order to obtain 
accurate lateral distribution functions. In the following 
we are going to analyze the influence of these parameters 
on the performance of the CE. 



The lower threshold E^^ should be chosen in the re- 
gion where the electromagnetic shower cannot be treated 
anymore as one-dimcnsional, and the lateral spread of 
electrons, positrons and photons becomes important. 
Fig. 5 shows the longitudinal and lateral profiles for 
different values of -EJ^Jf^. As of a threshold of 10 GeV, 
the profiles do not change significantly anymore. A very 
similar approach was used in reference [5] . There, a more 
complicated set of cascade equations was solved which in- 
volved also angular deviations from the showor-axis, and 
secondary particles below 10 GeV were followed in a MC 
method. Here we see that it is sufficient to treat the 
problem above 10 GeV in a purely longitudinal way. 



^ phad 

Analogously to the lower threshold in electromagnetic 
CE, the proper choice of the lower threshold -Ejjj^n de- 
pends on where the one-dimensional assumption is not 
valid anymore for the hadronic part of the shower. In 
order to test this we show in Fig. 6 the longitudinal 
and lateral and distribution function of muons, which 
are direct decay products of pions and kaons. A value 
of E^^ = W^GeV provides sufficient precision for both 
profiles. 



B. Tests on a single shower basis 

By evolving the high energy part of a shower with the 
MC-method, we are able to reproduce the natural fluc- 
tuations, which are primarily due to the varying depth 
of the first interactions. All particles which fall below 
a threshold fEo are used in the initial condition for the 
CE. In order to show that the CE are solved correctly 
for an arbitrary initial condition, we compare to the MC 
method by computing the high energy part in exactly 
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FIG. 5: Longitudinal profile of electrons and positrons (upper 
figure) and the lateral distribution function (lower figure) of 
electrons/positrons (lower curves) and photons (upper curves) 
for a photon induced shower and different values of E^-^^. 
In the lateral case, the four curves for > 3.16GeV are 

indistinguishable. 



the same way for both approaches. Technically this can 

be done by choosing the same seed for the pscudo ran- 
dom number generator in the computer program. Fig. 
7 shows such a comparison for the longitudinal and lat- 
eral profile of electrons and positrons. The thinning level 
for the MC method is 10~7. We see a slight sensitiv- 
ity of the longitudinal profile to the number of bins used 
in the numerical solution of the cascade equations. The 
shower maxima are at 738, 742, 739, and 740 g/cm? for 
the MC method, and the 10, 30 and 50 bin solution of 
the CE, respectively. The lateral distribution function is 
very insensitive to the number of bins. 



C. Statistical properties - fluctuations 

Next, we compare the statistical properties of two sets 
of proton induced lO^^eF showers calculated with CE 
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FIG. 6: The dependence of longitudinal and lateral profiles 

of muons as a function of -E^i^n. In the lateral case, the two 
curves for -E^JSn ^ W^GeV are indistinguishable. 




100 1000 
distance from axis [ml 



10000 



FIG. 7: An example of a shower calculated in the same way 

to down to 0.001-Eo, using below MC method (triangles) with 
10~^ thinning , or CE with different binnings. 



(500 showers) and MC-method (500 showers, 10~^ thin- 
ning). In Fig. 8 one sees the distribution of the shower 
maxima for the MC and CE methods. The two distribu- 
tions agree well. The threshold /, where CE takes the 
initial condition from the high energy MC part and com- 
putes the shower numerically, was chosen to be 0.001. 
The influence of the parameter / on the fluctuations is 

shown in Fig. 9, by plotting a = {X^^J - (X^ax)^ 
against /. The smaller the value of /, the further the 
initial portions of the shower is followed exactly rather 
than with the CE. We see that oven for / approaching 
unity, the fluctuations seen at small / are reproduced to 
a great extent. This shows that natural fluctuations arise 
for the most part from the depth of the flrst interaction 
of the cosmic ray in the atmosphere. 

D. Comparisons with CORSIKA 

Finally, we compare some of our results to CORSIKA 
simulations. CORSIKA [13] is a well tested simula- 
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FIG. 8: Comparison of a shower maximum distribution for 
10^^ eV proton induced showers. 



8 



70 



65 



E 60 



55 



50 



••. «i 



10'^ lO "* 10'^ 10"^ 10"^ 0.9 0.99 0.999 
f 



FI G. 9: The fluctu ations of the shower maximum a = 

\J (^max) — (^max)^ as a fuiiction of the fractional energy 
threshold. 1000 showers have been calculated for each value 
of/. 



tion package, which can be configured to use the same 
external models employed here for our hybrid model: 
QGSJET for the high energy hadronic part, GHEISHA 
for low energy hadrons and EGS4 for the electromagnetic 
part. Fig. 10 shows in the upper panel a comparison of 
the shower maximum for proton and iron induced 30- 
degree inclined showers. The shower proton curves agree 
nicely; the iron curve is slightly higher in our case. This 
might be due to differences in the computation of nucleus- 
nucleus cross sections, which are calculated from nucleon- 
nucleon cross sections using the Glauber method. We use 
./max ~ 0.001 < 1/56 which avoids calculating the hmc- 
tions Wmn for all 56 possible nuclei. The lower panel 
compares the average lateral distributions functions of 
5 X lO^^eV proton showers. They agree nicely for elec- 
trons/positrons and photons; compared to CORSIKA we 
produce slightly less muons, certainly due to the fact that 
we neglect photo-nuclear reactions at this stage. 



E. Summary 
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FIG. 10: Comparison to CORSIKA: the shower maximum as 
a function of the primary energy (top figure) for 30 degree 
inclined proton and Fe showers and the lateral distribution 
function for the average of ten 5 x 10^^ eV proton induced 
vertical showers (bottom figure). 



We introduced a hybrid approach to air shower simu- 
lations which uses cascade equations as proposed in ref- 
erence [8]. This method consists of applying traditional 
Monte Carlo methods where natural fluctuations or the 
lateral spread of particles are important. Particles are 
passed to the cascade equations via the initial condition. 
The low energy source function obtained from the cas- 
cade equations provides the probability distribution of 
low energy particles, whose further propagation is fol- 
lowed by Monte Carlo. The hybrid approach takes ad- 
vantage of the fast solutions of the cascade equations and 
provides detailed knowledge about each low energy par- 
ticle, such as position, energy and arrival time. 

Consistency checks have been made by comparing the 
hybrid CE approach to traditional Monte Carlo. The two 
methods agree nicely within a small error. The longitudi- 
nal profiles obtained with the CE approach are somewhat 
sensitive to the binning which enters in the numerical so- 
lution of the CE, if less than w 30 bins per decade are 
used. The lateral distribution functions are very stable 
against these technical parameters. 

The hybrid technique introduced here is faster than a 
traditional MC by at least a factor of 20 at lO^^eF. 
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